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A  percolation  model,  accounting  for  polydispersion  of  powders  and  presence  of  pore  formers  (i.e.  poros¬ 
ity),  is  presented  to  predict  effective  properties  of  composite  electrodes  for  solid  oxide  fuel  cells,  such 
as  the  three-phase  boundary  length  and  the  mean  hydraulic  radius.  Porosity  affects  both  numbers  of 
contacts  (so  probabilities  of  connection)  and  number  of  particles  per  unit  volume.  Both  these  effects, 
together  with  granulometric  distribution,  are  accounted  for  the  estimation  of  effective  properties.  As  a 
consequence,  the  theory  can  predict  numbers  of  contacts,  coordination  numbers  and  therefore  effective 
properties  of  the  electrode  for  multicomponent  polydisperse  mixtures. 

Model  simulations  show  that  the  three-phase  boundary  length  sharply  decreases  as  porosity  increases 
while  the  effects  of  polydispersion  of  powders  are  less  pronounced,  although  significant,  suggesting  that 
these  features  should  be  considered  in  SOFC  electrode  models. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Porous  composite  electrodes  are  used  in  the  architecture  of 
solid  oxide  fuel  cells  (SOFCs)  [1,2].  A  composite  electrode  is  a 
sintered  random  packing  of  electron-conducting  particles  and  ion¬ 
conducting  particles  where  electrochemical  reactions  of  oxidation 
or  reduction  occur.  As  an  example,  a  typical  composite  cathode 
is  a  mixture  of  lanthanum-doped  strontium  manganite  (LSM)  as 
electronic  conductor  and  yttria-stabilized  zirconia  (YSZ)  as  ionic 
conductor  [3].  The  charge  transfer  step  is  considered  to  take  place 
where  the  reaction  constituents  (i.e.  electrons,  ions  and  chemical 
species)  can  coexist.  These  sites  are  the  three-phase  boundaries 
(TPBs)  where  the  distinct  phases,  transporting  the  three  participat¬ 
ing  species,  meet. 

The  electrode  performance  is  strictly  dependent  on  its  morphol¬ 
ogy.  In  particular  the  extension  of  the  reactive  zone,  the  effective 
conductivities  of  the  phases  and  the  morphological  properties 
related  to  the  gas  phase  are  function  of  particle  size,  composition, 
porosity  and  sintering  conditions.  As  a  consequence,  the  engineered 
design  of  electrodes  requires  understanding  and  optimization  of 
their  morphology  [3,4]. 

The  main  aim  of  studying  the  electrode  morphology  is  to  predict 
the  effective  functional  properties  (e.g.  the  mean  pore  diameter, 


*  Corresponding  author.  Tel.:  +39  3334687172;  fax:  +39  (0)502217866. 

E-mail  addresses:  antonio.bertei@for.unipi.it,  berteiantonio@libero.it  (A.  Bertei). 

0378-7753 /$  -  see  front  matter  ©  2011  Elsevier  B.V.  All  rights  reserved, 
doi:  10.1 01 6/j.jpowsour.201 1.06.087 


the  effective  conductivity  of  phases,  the  three-phase  boundary 
length)  starting  from  measurable  parameters  such  as  the  granulo¬ 
metric  distribution  of  the  powders  and  their  composition.  Methods 
for  characterizing  electrode  morphology  can  be  broadly  classi¬ 
fied  as:  (i)  experimental  techniques,  (ii)  computer  simulations 
of  random  packings  and  (iii)  predictions  based  on  percolation 
theory. 

Microstructural  parameters  can  be  measured  on  existing  sam¬ 
ples  by  image  analysis  on  scanning  electron  microscopy  (SEM) 
images  or  by  other  specific  experimental  techniques  (e.g.  measure¬ 
ments  of  effective  conductivity  or  BET  adsorption  measurements). 
The  morphology  of  the  packing  can  be  reconstructed  (e.g.  by  ion 
beam  SEM)  and  analyzed  [5-7].  These  techniques  allow  correla¬ 
tions  of  electrode  properties  to  powder  characteristics,  but  they 
rely  on  the  resolution  of  the  image  and/or  phase  identification, 
making  analysis  of  small  features  difficult  and,  in  some  cases,  uncer¬ 
tain. 

Composite  electrodes  are  often  modeled  as  random  mixtures  of 
spherical  particles  containing  at  least  two  solid  phases:  electron¬ 
conducting  phase  and  ion-conducting  phase  [8-10].  The  electrode 
structure  can  be  reproduced  by  using  a  packing  algorithm  and 
morphological  parameters  are  evaluated  on  the  simulated  packing 
[11-13].  Despite  crude  simplifications,  this  method  is  faster  and 
may  be  more  accurate  than  experimental  techniques  since  it  is  not 
affected  by  experimental  and  image-processing  uncertainties. 

Microstructural  properties  can  also  be  evaluated  by  using 
a  stochastic  theory  regarding  the  particle  arrangement  in  the 
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packing,  i.e.  by  using  percolation  theory.  This  mathematical  tool 
has  been  widely  used  to  describe  morphological  properties  of  ran¬ 
dom  packings  of  solid  particles,  in  particular  in  SOFC  modeling 
[8,9,14],  due  to  its  easier  applicability.  Unlike  simulation  methods 
based  on  numerical  packing  generation,  it  is  a  predictive  method, 
since  functional  properties  are  obtained  directly  from  measurable 
starting  input  parameters,  such  as  particle  size  and  powder  com¬ 
position. 

Percolation  theory  relies  on  the  prediction  of  coordination  num¬ 
bers  of  particles,  i.e.  the  number  of  contacts  made  by  a  particle 
with  its  neighboring  particles.  Functional  morphological  properties 
are  predicted  starting  from  the  evaluation  of  coordination  num¬ 
bers.  Several  percolation  models  have  been  developed  [14-16], 
based  on  different  basic  relationships  used  to  calculate  coordi¬ 
nation  numbers  as  a  function  of  particle  size  and  composition. 
In  a  previous  paper  [17]  the  authors  showed  a  critical  compari¬ 
son  of  existing  percolation  theories  and  proposed  a  new  method 
based  on  the  improvement  of  the  Suzuki  and  Oshima  theory  [15] 
to  obtain  a  better  agreement  with  simulated  and  experimen¬ 
tal  results  and  an  extension  of  the  theory  to  multicomponent 
polydisperse  powders.  As  the  basis  for  the  development  pre¬ 
sented  in  this  work,  the  essential  features  of  the  new  method  are 
recalled  in  the  background  of  the  theoretical  section  (see  Section 
2.1). 

The  present  paper  develops  an  extension  of  the  previous  the¬ 
ory  to  systems  where  porosity  is  controlled  by  addition  of  pore 
formers  in  the  particle  mixture.  Pore  formers  are  particles  of 
materials  which  decompose  at  high  temperatures,  mixed  with 
the  powders  before  sintering  in  order  to  increase  the  porosity 
of  the  packing  after  sintering.  Typical  pore  former  materials  are 
starch  [18,19]  and  graphite  [20,21]  which  burn  in  oxidizing  atmo¬ 
sphere  during  the  sintering.  Although  commonly  used  in  electrodes 
fabrication,  pore  formers  have  not  considered  yet  in  percolation 
theories. 

The  focus  of  this  paper  is  then  on  the  application  of  the  new 
percolation  model  to  predict  effective  properties  for  composite 
electrodes  of  SOFCs,  showing  in  particular  the  effects  of  both  the 
presence  of  pore  formers  and  the  polydispersion  of  powders.  In  the 
theoretical  section,  the  percolation  model  is  first  presented  (Sec¬ 
tion  2.1)  and  then  the  extension  to  consider  pore  former  effects 
is  developed  (Section  2.2).  After  a  section  dedicated  to  analyze  if 
overlapping  effects  affect  the  percolation  theory  estimations  (Sec¬ 
tion  2.3),  the  derivation  of  effective  functional  properties  based 
on  percolation  theory  is  shown  (Section  2.4).  Model  outcomes  are 
specifically  discussed  in  chapter  3. 


2.  Theory 

2.1.  Background 

This  section  summarizes  the  key  features  of  the  percolation 
model  previously  developed  [17].  This  model  was  compared  and 
validated  with  other  percolation  models  [14-16]  over  a  wide  range 
of  particle  sizes  and  compositions.  Despite  in  electrode  fabrication 
only  two  conducting  phases  are  used  (i.e.  electron-conducting  and 
ion-conducting),  also  a  third  solid  phase  is  here  considered,  which 
may  represent  pore  former  particles  (as  discussed  in  Section  2.2) 
or  other  types  of  particles. 

Let  us  consider  a  ternary  polydisperse  mixture  of  rigid  spher¬ 
ical  particles,  namely  k,  1  and  f.  The  packing  is  characterized  by  a 
distribution  of  mk  particle  sizes  for  k-particles  (i.e.  different  radii 
rk\y  rk2. .  .rkm  ),  m/  for  1-particles  and  m^for  f-particles.  For  a  generic 
kind  of  hi-particles  (h  =  k,  1  or  f  while  i  varies  from  1  to  mk,  mi  or 
mf),  characterized  by  a  radius  rhi,  basic  definitions  of  ratio  of  radii 
with  another  generic  kind  of  gj-particles  (g  =  k,  1  or  f  while  j  varies 


from  1  to  mk ,  m/  or  mf),  numerical,  volumetric  and  surface-specific 
fractions  before  sintering  are  as  follows: 


Phi,gj  — 


rhi 


(1) 


Zhi  = 


x^mk ; 


't'hi/rl 


ETi  *»/*%  +  E;Ji  i'lj/rfj  +  Ej7i  Vrj 


(2) 


Von  = 


£hirhi 


(3) 


Ski  — 


%hirM 


Ej=i  Skjrlj  +  Ei= 1 4ry  +  Ej=fi  &rj 

= _ tyhih I, _ 

Eh i  ^kj/rkj  +  Ej=i  VV  n +  Ej7i  ^s/rs 


(4) 


where  nvhi  represents  the  number  of  hi-particles  per  unit  volume. 

The  packing  homogeneity  depends  on  the  values  of  Phigj.  For 
a  couple  of  particles  hi-gj,  segregation  occurs  if  the  smaller  par¬ 
ticles,  for  example  hi-particles,  can  pass  through  the  void  left  by 
three  larger  particles,  for  example  gj-particles,  in  contact  with  each 
other  [22].  Indeed,  by  joining  the  centers  of  bigger  gj-particles 
to  form  an  equilateral  triangle,  by  geometrical  considerations  the 
void  left  in  the  middle  can  contain  at  the  maximum  a  particle  of 
radius  (2/V3  -  1  )rs-:  if  rhi  is  higher,  segregation  cannot  occur.  This 
consideration  could  be  repeated  if  hi-particles  are  bigger  than  gj- 
particles,  providing  a  restriction  on  rs  as  rgj  >  (2/V3  -  1  )rhi.  Thus, 
a  sufficient  condition  to  avoid  segregation  for  the  multicompo¬ 
nent  polydisperse  mixture  is  given  if  Phigj  is  in  the  range  from 

pmin  =  2/V3  -  1  toPmax  =  (2/V3  -  l)-1  for  each  couple  of  particles 
hi-gj. 

The  generic  average  number  of  contacts  between  hi-  and  gj- 
particles  Zhigj,  i.e.  the  number  of  gj-particles  in  contact  with  a  hi- 
particle,  is  proposed  to  be  proportional  to  surface-specific  fraction 
of  gj-particles  as  in  Suzuki  and  Oshima  model  [15]: 

Zhi,gj  —  Sgj  Zhi,gj\r  ( 5 ) 


Zhi  gj  represents  the  limit  of  the  number  of  contacts  that  a 

hi-particle  makes  with  gj-particles  in  a  binary  hi-gj  mixture  when 
there  are  only  a  few  hi-particles  in  the  packing.  It  depends  only  on 
the  ratio  of  radii  Phigj  and  it  is  estimated  as  follows: 


'  3(2  -  \/3)(P)„- ^  +  1 ) 

0+Phi,gi-(Phi,s(Phi,s  +  2)f5} 
3(2-^)(p-1g).  +  l) 


if  Pw,a  >  1 

i^Phi,gj  <  1 


(6) 


Eq.  (6)  represents  the  improvement  of  Suzuki  and  Oshima  theory 
as  extensively  discussed  in  [17]. 

For  a  generic  hi-particle,  the  number  of  contacts  between  hi  and 
all  kinds  of  g-particles  Zhig,  where  g  represents  k-,  1-  or  f-particles, 
is  calculated  as  the  summation  of  the  number  of  contacts  that  a 
hi-particle  makes  with  other  gj-particles: 

mg 

zhi, g  =  Yfa  withg  =  k,  lorf  (7) 

j= i 
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Fig.  1.  Illustration  of  a  composite  electrode  with  A-(i.e.  entirely  connected),  B-(i.e. 
partly  connected)  and  C-(i.e.  isolated)  clusters  indicated  for  k-phase. 


The  average  coordination  number  of  hi-particles  (Zw)  is  then 
defined  as  the  sum  of  the  average  number  of  contacts  that  hi- 
particles  make  with  k-,  1-  and  f-particles: 


Zhi  —  Zhi,k  +  Zhi,l  +  Z hij  (8) 

The  overall  average  number  of  contacts  for  all  h-particles  Zh>h, 
i.e.  the  average  number  of  contacts  that  a  h-particle  makes  with 
other  homologue  h-particles,  is  proposed  to  be  equal  to  the  total 
number  of  h-h  contacts  in  the  mixture  divided  by  the  total  number 
of  h-particles  as  already  proposed  by  Chen  et  al.  [14]: 

v  ^-1  ShiZhi.h  .  ,  ,  i  i  f 

Zh,h  =  ^mh  -  withh  =  k,  lorf  (9) 

2^1=1  Shi 

As  shown  in  Fig.  1 ,  a  percolating  (or  A-)  cluster  is  a  cluster  of  par¬ 
ticles  of  the  same  type  that  extends  through  the  entire  thickness  of 
the  electrode.  On  the  other  hand,  B-  and  C-clusters  are  not  entirely 
connected  because  B-clusters  are  short  networks  connected  only  to 
the  respective  collector  while  C-clusters  are  completely  isolated  [9]. 
The  probability  of  connection  ph ,  i.e.  the  probability  that  a  generic 
h-particle  belongs  to  a  percolating  cluster  of  the  same  kind  of  par¬ 
ticles,  depends  only  on  the  overall  average  number  of  contacts 
among  h-particles  Zh  h  [16].  The  probability  of  connection  ph  is  then 
calculated  as  [17]: 


Ph  =  1  - 


/  4.236  —  Zhh 
^  2.472 


3.7 


(10) 


Eq.  (10)  can  also  be  used  to  individuate  percolation  thresholds, 
i.e.  the  minimum  volume  fraction  of  particles  which  determines 
percolating  clusters  of  particles  of  the  same  type.  In  particular,  tak¬ 
ing  generically  h-particles  as  reference,  when  ph  equals  0  (i.e.  for 
Zh  h  =  1 .764)  in  the  mixture  there  are  no  more  percolating  clusters  of 
h-particles.  It  means  that  the  h-phase  is  no  more  connected  through 
the  whole  thickness  of  the  packing  and,  as  a  result,  the  effective 
conductivity  of  this  phase  approaches  0. 


2.2.  Extension  of  the  percolation  model  including  pore  formers 
effects 


The  percolation  theory  presented  in  the  previous  section  is 
extended  to  take  into  account  the  presence  of  pore  formers  in  the 
particle  mixture.  Pore  former  particles  are  not  present  in  the  final 
structure  because  they  decompose  during  the  sintering. 


In  a  ternary  mixture  of  k-,  1-  and  f-particles  as  considered  in 
Section  2.1,  where  f-particles  represent  pore  former  particles,  we 
assume  that  after  sintering  f-particles  disappear  leaving  holes  of 
the  same  shape  and  dimension  of  the  particle  of  pore  former  that 
occupied  that  position  before  sintering.  The  resulting  structure, 
obtained  after  sintering,  will  contain  only  k-  and  1-particles.  It  is 
then  assumed  that  size  and  position  of  each  k-  and  1-particle  do 
not  change  during  the  sintering.  This  assumption  is  valid  only  if  the 
structure  does  not  collapse  during  the  sintering.  Structure  collapse 
may  reasonably  occur  for  high  volume  fractions  of  pore  formers 
(such  as  to  create  porosity  higher  than  50-60%)  or  if  pore  former 
particles  are  much  bigger  than  k-  and  1-particles,  resulting  in  big 
holes  whose  surrounding  structure  is  mechanically  unstable.  Even 
when  pore  former  particles  are  much  smaller  than  k-  and  1-particles 
the  structure  may  collapse  since  the  smaller  f-particles  surround 
other  coarser  particles  isolating  them  [23  ] :  during  the  sintering,  the 
gradual  disappearance  of  pore  former  particles  brings  the  particles 
placed  above  them  to  come  down. 

In  order  to  calculate  numbers  of  contacts,  Eqs.  (l)-(4)  remain 
valid  for  all  the  types  of  particles  (i.e.  k,  1  and  f)  because  they 
represent  the  ratio  or  radii  and  numerical,  volumetric  and  surface- 
specific  fractions  before  sintering.  These  quantities  still  continue 
to  have  a  physical  meaning  after  sintering  since  f-particles  (and 
their  correlated  quantities)  represent  the  spherical  holes  left  in  the 
structure.  Nevertheless,  after  sintering  Eqs.  (5)-(7)  and  (9)-(10) 
continue  to  make  sense  only  for  k-  and  1-particles  (i.e.  for  h  and 
g  equal  to  k  or  1),  despite  numerical  and  surface-specific  fractions 
involved  in  calculations  must  be  determined  by  using  Eqs.  (l)-(4) 
where  pore  formers  are  considered.  Numbers  of  contacts  regard¬ 
ing  f-particles  do  not  have  any  physical  meaning  after  sintering, 
although  f-particles  affect  the  determination  of  numbers  of  con¬ 
tacts  of  k-  and  1-particles  since  they  represent  missing  particles. 
Coordination  numbers  for  k-  and  1-particles  can  be  calculated  by 
using  Eq.  (8),  nevertheless  the  contacts  with  f-particles  must  be 
ignored: 

zki  —  zki,  k  +  Zfa  i  (11) 

Zli  =  Zli,k  +  Zli,l  (12) 

Pore  formers  affect  not  only  the  calculation  of  numbers  of  con¬ 
tacts  but  also  the  number  of  particles  per  unit  volume.  The  final 
porosity  c//™  of  the  structure  (i.e.  after  sintering)  can  be  calculated 
as: 

0f  =  </>">  + 4  (13) 

where  </^n  is  the  initial  porosity  of  the  mixture  (i.e.  before  sinter¬ 
ing)  when  all  k-,  1-  and  f-particles  are  present.  c//g  is  the  fraction  of 
porosity  added  by  pore  formers  calculated  as: 

4  =  nv’bs^2^fi{  4/3)7tr^  (14) 

1=1 

where  nv’bs  is  the  number  of  particles  (k,  1  and  f)  per  unit  volume 
(voids  included)  before  sintering  estimated  as: 

nv,bs\  _  _ 1  ~  _ 

Eh  4Tfcf(4/3)7rr3  +  ^  ^(4/3)^  +  ^(4/3>rrj 

(15) 


The  number  of  particles  per  unit  volume  after  sintering  nv>as  is 
calculated  as: 


^v,as  _  ftV,bs 


mk 


mi 


+ X/' 


(16) 


9432 


A  Bertei,  C.  Nicolella  /  Journal  of  Power  Sources  196(2011)  9429-9436 


The  initial  porosity  </^n  depends  on  composition  and  gran¬ 
ulometric  distribution  of  powders  before  sintering  but  this 
relationship  is  not  well  known  for  multicomponent  mixtures.  Thus, 
if  no  particular  measurement  is  available,  a  very  rough  estimation 
is  =  0.36,  i.e.  the  porosity  of  a  random  close  packing  of  mono¬ 
sized  spheres  of  the  same  type  [16,24-26].  In  general,  for  binary 
mixtures  of  monosized  (i.e.  not  polydisperse)  particles  with  differ¬ 
ent  radii,  initial  porosity  can  be  evaluated  with  more  accuracy  as 
reported  by  Ben  Aim  and  Le  Goff  [27]. 


2.3.  Overlapping  of  particles 

In  percolation  theory,  particles  are  considered  as  spheres  that 
do  not  overlap  as  they  are  rigid.  On  the  other  hand,  in  the  man¬ 
ufactory  of  electrodes  for  SOFCs,  particles  are  compressed  and 
sintered.  Then,  assuming  for  instance  that  before  sintering  parti¬ 
cles  are  spherical,  after  that  operation  we  should  expect  that  their 
shape  shows  overlapping  regions.  This  implies  that  some  properties 
of  the  packing,  in  particular  numbers  of  contacts  and  coordination 
numbers  of  particles  but  also  the  number  of  particles  per  unit  vol¬ 
ume,  could  differ  from  the  estimation  made  by  assuming  particles 
as  rigid  spheres. 

The  discrepancy  between  the  properties  of  the  sintered  pack¬ 
ing  (i.e.  with  overlapping)  and  the  modeled  packing  as  considered 
by  percolation  theory  (i.e.  without  overlapping)  increases  when 
contact  angles  among  particles  increase.  The  graphical  definition 
of  contact  angles  between  a  k-  and  a  1-particle  (i.e.  0ki  and  0ik)  is 
reported  in  Fig.  2. 

In  order  to  estimate  the  effects  of  overlapping  at  least  on  the 
number  of  particles  per  unit  volume,  let  us  consider  a  binary  ran¬ 
domly  packed  mixture  of  k-  and  1-particles  not  polydisperse.  As 
already  mentioned  by  Janardhanan  et  al.  [28],  we  can  assume  that: 


1.  Particles  are  spherical  before  sintering 

2.  The  sintering  process  creates  overlapping  regions  among  parti¬ 
cles  characterized  by  mean  contact  angles  0kk ,  0kl,  0ik  and  Op 

3.  After  sintering  particles  keep  quite  the  same  position,  same 
radii  rk  and  r/  and  the  spherical  shape  except  in  the  overlap¬ 
ping  region:  the  contact  region  between  two  particles  is  a  perfect 
circle  of  radius  rk  sm0ki  =  ri  sm0ik 

4.  Presence  of  overlapping  regions  made  by  three  or  more  particles 
is  neglected 


We  can  calculate  the  fraction  of  volume  lost  in  each  kind  of  contact. 
We  call  v!ost,ij  the  volume  of  a  generic  i-particle  lost  in  a  contact  i-j 
(characterized  by  an  angle  Oij).  According  to  Fig.  2,  the  lost  volume  is 


equal  to  the  volume  of  a  spherical  segment  of  one  base  with  height 
r,-(  1  -  cos  Gy),  so  we  obtain: 

'  KSt’kl  =  JrkO  ~  cos0kif  [3  rk  -  rk(  1  -  cos#w)] 

yhstM  =  _  cos0k/(f  [3rfc  _  rk(  i  _  cos#tt)] 

viost,ik  =  |r2(1  _  cos 6lkf  [3 r,  -  r,(  1  -  cos#,k)] 

k  v/0*"  =  |r2(  1  -  cos#,,)2  [3r,  -  r,(l  -  cos#,,)] 

with  the  constraint  (see  point  3  mentioned  above): 

rksin  0ki  =  risin0ik  (18) 


that  links  together  0ki  and  0[k. 

Assuming  that  percolation  properties  of  the  mixture  can  be 
estimated  by  using  percolation  theory  described  in  Section  2.1 
although  particles  are  not  rigid  spheres,  the  total  volume  lost  in 
the  contacts  for  each  particle,  k  or  1,  is: 

VkSt  =  zuK°St’kl  +zk,kV'kSt'kk  (19) 


i /lost  7  I/I 
vl  =zl,kvi 


lost,  Ik 


+  zuvl 


lost, ll 
l 


(20) 


Thus,  the  effective  volumes  of  the  particles,  V(f  and  V'!n ,  are: 


,'ff 


Vf  =  (4/3)m3k  -  V i 


lost 


yf  =  (4/3  )jnf  -  Vj 


lost 


(21) 

,  -  v,  (22) 

These  relationships  are  useful  to  calculate  the  total  number  of 
particles  per  unit  volume  (voids  included)  with  overlapping  regions 

jjV,  overlap . 

i-# 


^v, overlap  _ 


ftvf  +  (i-4)y,e// 


(23) 


If  Oy  >  0,  nv’overlap  will  be  higher  than  the  value  of  nv’bs  that  we 
could  calculate  without  considering  overlapping  effects  as  in  Eq. 
(15). 

For  reasonable  contact  angles,  i.e.  0kk  =  Op  =  ma x(0ki,0ik)  <  1 5°  as 
commonly  accepted  in  modeling  works  [8,9,14],  the  errors  that  we 
make  on  nv’overlap  if  we  do  not  consider  the  overlapping  effects  (i.e. 
using  (4/3)7rr3  instead  of  Vf^)  will  be  less  than  0.52%  (the  worst 
situation  is  represented  by  a  binary  mixture  of  k-  and  1-particles 
of  the  same  size).  It  means  that  overlapping  effects  are  negligible 
in  the  calculation  of  nv’overlap‘,  secondarily,  it  could  mean  that  also 
the  approximation  of  the  mixture  as  a  packing  of  rigid  spheres  is 
reasonable,  so  we  can  estimate  coordination  numbers  of  particles 
and  other  properties  by  using  the  relationships  of  percolation  the¬ 
ory  presented  in  Sections  2.1  and  2.2  even  if  particles  are  partly 
overlapped. 


2.4.  Calculation  of  functional  parameters 


As  far  as  modeling  of  electrodes  for  SOFCs  is  concerned,  the  most 
relevant  morphological  parameters  are  the  three-phase  bound¬ 
ary  (TPB)  length  per  unit  volume  XPPB,  the  contact  area  between 
connected  k-  and  1-particles  per  unit  volume  avkl,  the  surface  area 
exposed  by  k-  (or  1-)  particles  to  gas  phase  per  unit  volume  avk  (or 
a[)  and  the  mean  hydraulic  radius  Rh  for  the  description  of  trans¬ 
port  of  gaseous  products  or  reagents.  The  method  proposed  for  the 
estimation  of  these  parameters  is  presented  in  the  following. 


2.4 A.  TPB  length 

When  considering  two  types  of  conducting  particles,  k  and  1  (e.g. 
k,  electron-conducting  phase,  1,  ion-conducting  phase),  the  TPB  is 
the  contact  perimeter  between  overlapping  k-  and  1-particles.  The 
third  phase  is  the  gas  phase,  generally  assumed  entirely  connected 
(i.e.  we  can  connect  two  generic  points  of  the  gas  phase  passing 
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only  within  the  gas  phase)  for  porosity  higher  than  25%  [  1 1  ].  At  the 
TPB,  electrons,  ions  and  chemical  species  can  simultaneously  coex¬ 
ist.  A  k-1  contact  will  be  active  for  the  charge-transfer  reaction  only 
if  the  k-particle  belongs  to  a  k-path  connected  to  the  k-collector 
and  the  1-particle  to  a  1-path  connected  to  the  1-collector.  Only  in 
this  situation  reaction  can  occur  at  the  TPB,  e.g.  in  a  cathode  elec¬ 
trons  flow  from  the  electrode  interconnect  to  the  TPB  to  react  with 
gaseous  oxygen  while  oxygen  ions  flow  from  the  reaction  site  at  TPB 
to  the  dense  electrolyte.  The  contributions  of  B-clusters  (i.e.  con¬ 
tacts  between  A-B,  B-A  or  B-B  clusters  of  different  kinds,  see  Fig.  1 ) 
can  generally  be  neglected,  because,  within  percolation  thresholds, 
their  length  is  very  low  when  the  thickness  of  the  electrode  is  rea¬ 
sonably  big  if  compared  with  mean  particle  size  [9].  In  this  case, 
only  the  contacts  among  particles  belonging  to  A-clusters  are  taken 
into  account. 

The  TPB  length  per  unit  volume  is  calculated  starting  from  the 
TPB  length  per  ki-lj  contact  27nT,  i.e.  the  contact  perimeter,  where  rf. 
is  the  neck  radius  of  the  circle  formed  by  the  ki-lj  contact,  calculated 
as: 

r-j  =  min(rfc(,  rtj)  sin(6>?)  (24) 

coupled  with  the  definition  of  the  contact  angle  0?.: 

=  max(0w>!/,  fyirf)  (25) 

The  neck  radius  is  calculated  by  referring  to  the  smaller  par¬ 
ticle  between  ki  and  lj,  i.e.  by  using  the  minimum  particle  radius 
between  rki  and  ry  and  the  maximum  contact  angle  between  Okiij 
and  6ij>ki. 

The  TPB  length  per  unit  volume  is  calculated  through  the  fol¬ 
lowing  procedure: 

1.  The  contact  perimeter  2tt rc.  is  multiplied  by  the  number  of  ki-lj 

contacts  per  unit  volume,  equal  to  nv’bs^kiZki  ij 

2.  The  sum  is  extended  to  the  mk  mi  couples  of  contacts  between 
k-  and  1-particles 

3.  Only  particles  belonging  to  connected  clusters  are  considered  by 
multiplying  by  pk  and  pi 

Accordingly,  the  TPB  length  per  unit  volume  is  calculated  as: 


contact  perimeter  27rf.  for  each  contact  ki-lj,  we  take  into  account 
the  contact  area  7r(rf.)2: 


mk  (  ml 

ah  =  ""’WE  I 

>=i  \  J=i 

This  last  parameter  characterizes  reaction  sites  if  the  reaction 
takes  place  in  the  whole  surface  of  contact  between  k-  and  1- 
particles  instead  of  the  contact  perimeter. 

2.4.3.  Surface  area 

The  total  surface  area  exposed  by  k-particles  to  gas  phase  per 
unit  volume  avk  has  a  role  concerning  adsorption  (or  desorption) 
of  substances  from  gas  phase  to  k-particles.  It  is  calculated  by 
multiplying  the  total  number  of  particles  per  unit  volume  before 
sintering  nv’hs  with  the  surface  area  of  a  ki-particle  4 7zrfd  and 
extending  the  summation  to  the  mk  particle  sizes: 

mk 

al  =  47rn'’-te^4  (28) 

1=1 

Eq.  (28)  can  be  applied  also  for  1-particles.  The  expression  takes 
into  account  the  total  external  surface  of  k-particles  without  con¬ 
sidering  the  fraction  of  surface  lost  in  contacts.  Considering  contact 
angles  of  15°,  due  to  the  fraction  of  surface  lost  in  contacts,  Eq.  (28) 
overestimates  avk  in  the  order  of  1 1 .4%  in  the  worst  case  of  a  binary 
mixture  of  monosized  k-  and  1-particles  with  rk  =  ri  without  pore 
formers. 

2.4.4.  Hydraulic  radius 

Referring  to  the  gas  phase,  the  main  parameter  is  the  mean 
hydraulic  radius  Rh,  estimated  according  to  Bird  et  al.  [29]  as  the 
ratio  of  the  volume  available  for  gas  passage  per  unit  volume,  i.e. 
the  final  porosity  </^n,  and  the  total  solid  surface  exposed  to  gas 
phase  per  unit  volume,  equal  to  the  sum  of  avk  and  av{  both  calcu¬ 
lated  according  to  Eq.  (28).  By  this  definition,  the  hydraulic  radius 
is  calculated  as: 


Rh  = 


4n 


4  (££&»£ +  £M) 


(29) 


mk  /  ml 

XTPB  =  nv'bsPkPi2xJ2  (  tki'ffZkijjry 

1=1  \  j=l 

In  particular,  in  Eq.  (26)  are  the  numerical  fractions  before 
sintering  (i.e.  considering  also  the  presence  of  pore  former  par¬ 
ticles)  and  nv>hs,  Zki y,  pk  and  pi  are  calculated  according  to  the 
relationships  presented  in  Sections  2.1  and  2.2. 

All  the  contact  angles  0f.  (i.e.  mk  mi  parameters)  must  be 
assigned.  It  is  reasonable  to  assume  a  single  value  for  all  of  them 
and  call  it  0:  in  this  situation  Eq.  (26)  is  similar  to  the  expression 
already  proposed  by  Chen  et  al.  [14].  According  to  the  considera¬ 
tions  formulated  above  in  Section  2.3,  0  must  be  assumed  in  the 
order  of  15°  or  less  in  order  to  apply  the  relationships  of  percola¬ 
tion  theory.  Other  values  of  0  could  be  used,  referred  to  different 
sintering  conditions,  but  for  big  contact  angles  (i.e.  0  >  30° )  particles 
are  no  more  spherical,  so  percolation  theory  cannot  be  used. 


2.4.2.  Contact  area 

In  a  similar  way,  we  can  evaluate  avkV  i.e.  the  contact  area 
between  connected  k-  and  1-particles  per  unit  volume.  The  calcula¬ 
tion  of  avkl  is  similar  to  that  of  X”pB,  only,  instead  of  considering  the 


Since  the  fractions  of  surface  lost  in  contacts  are  not  considered 
in  Eq.  (29),  as  previously  explained  for  ak,  Rh  is  underestimated  in 
the  order  of  1 0.2%  as  maximum  for  contact  angles  of  1 5°.  By  defini¬ 
tion,  the  mean  diameter  of  pores  is  equal  to  4  times  the  hydraulic 
radius. 


3.  Results  and  discussion 


The  model  outcomes  are  the  functional  parameters  XPPB,  avkp  avk 
and  Rh,  calculated  as  described  in  the  theoretical  section  (Section 
2.4).  In  order  to  illustrate  the  model  potentiality,  functional  param¬ 
eters  are  calculated  for  several  case-studies  in  different  conditions 
to  explain  the  role  of  starting  input  parameters  (i.e.  radii  of  par¬ 
ticles,  composition,  porosity,  etc.)  on  final  morphological  results. 
In  particular,  special  attention  is  taken  in  illustrating  the  effects 
of  porosity  and  of  polydispersion  of  particles  as  the  new  features 
introduced  by  the  model.  All  results  are  referred  to  0  =  15°,  i.e.,  in 
an  accurate  notation,  0  =  0kk  =  0\\  =  =  15°.  In  the  following,  for 

numerical  and  volumetric  compositions  after  sintering  we  mean: 


i-as  _ 

Shi  ~ 


mk  ml 

and  $  =  ^  = 


(30) 
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Volume  fraction  ofk-particles,^°s 

Fig.  3.  TPB  length  per  unit  volume  as  a  function  of  volumetric  fraction  of  k-particles 
at  different  radii  of  particles  for  a  binary  not  polydisperse  mixture  without  pore 
formers  (0jn  =  0.36,  9  =  15°). 


mk  mi 

V,Hl  Jr  ■  and^k  =  ^  =  Yj**  (31) 

Z^I=1  ^  1  =  1  i=l 

for  k-  and  1-phases  (i.e.  h  =  k  or  1)  where  fa  and  i frhi  are,  respectively 
numerical  and  volumetric  fractions  before  sintering  as  already 
defined  in  the  theoretical  section  (Section  2.1). 

3.1.  Particle  size 

Let  us  consider  a  binary  mixture  of  monosized  (not  polydis¬ 
perse)  particles  without  pore  formers;  we  assume  for  simplicity 
a  constant  final  porosity  of  36%  for  each  composition  and  particle 
size  ratio.  Fig.  3  shows  XVTPB  as  a  function  of  t/^s  for  several  radii 
of  particles  for  0  =  15°.  When  particles  have  the  same  radius,  XVTPB 
reaches  a  maximum  at  t/^s  =  0.5  and  approaches  0  at  the  perco¬ 
lation  thresholds.  For  a  given  composition,  bigger  particles  yield 
lower  values  of  XvTpB  than  smaller  particles.  If  1-particles  are  bigger 
than  k-particles,  the  compositions  where  XVTPB  approaches  0  and  at 
which  XjPB  is  maximum  move  towards  lower  since  bigger  par¬ 
ticles  (1)  percolate  worse  than  smaller  particles  (k)  [23].  Note  also 
that  the  width  of  the  range  of  composition  at  which  X”pB  is  not  zero 
decreases  as  the  ratio  of  radii  Pkt  is  not  1. 

3.2.  Porosity 

Fig.  4  shows  in  solid  lines  XVJPB  as  a  function  of  t/^s  for  differ¬ 
ent  porosities  for  a  ternary  not  polydisperse  mixture  of  particles 
with  rk  =  ri  =  rf=0A  pun,  where  f-particles  represent  pore  former 
particles;  0  =  15°  is  used.  When  porosity  increases,  X^pB  decreases 
and  the  range  of  compositions  corresponding  to  XvTpB  >  0  becomes 
smaller. 

The  main  feature  of  our  percolation  model  is  that  the  presence 
of  pore  formers  has  a  double  effect  on  XVTPB:  first  it  enters  in  the  cal¬ 
culation  of  the  numbers  of  contacts  and  connection  probabilities, 
then  it  affects  the  calculation  of  number  of  particles  per  unit  volume 
after  sintering  (Eq.  (16)).  In  other  percolation  models  [9,14],  whose 
results  are  represented  in  Fig.  4  as  dotted  lines,  the  first  effect  is 
ignored,  yielding  an  overestimation  of  TPB  length  and  of  the  range 
where  XVTPB  >  0  when  pore  formers  are  used.  Provided  that  it  is  rea¬ 
sonable  that  pore  formers  reduce  the  number  of  contacts  between 
k-  and  1-particles,  our  percolation  model  appears  more  reliable  than 
other  presented  in  Fig.  4. 

Porosity  has  also  effects  on  the  mean  hydraulic  radius  Rh.  Fig.  5 
shows  Rk  as  a  function  of  composition  for  different  porosities  with 
rf=0A  [Jim,  rk  =  0A  pan  and  q  =  0.15  pm  and  with  0=15°.  Initial 
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Fig.  4.  TPB  length  per  unit  volume  vs  composition  at  different  porosities  for  mono¬ 
sized  particles:  comparison  between  our  model  (solid  lines)  and  other  models 
(dotted  lines)  [9,14],  in  which  pore  formers  do  not  affect  numbers  of  contacts 
(r/<  =  r/  =  r/=0.1  fxm,  6=  15°). 


porosity  is  considered  constant  and  equal  to  0.36  for  each  case. 
When  porosity  or  composition  of  the  bigger  particles  (1-  in  this  case) 
increase,  also  Rh  increases. 

3.3.  Pore  formers  size 

The  dimension  of  pore  former  particles,  i.e.  rp  is  another  impor¬ 
tant  parameter.  Given  the  final  porosity  equal  to  0.45  (i.e.  fixed 
i/ff=  0.140625)  with  initial  porosity  equal  to  0.36  and  kept  con¬ 
stant,  the  radius  rj  is  varied  in  a  ternary  mixture  of  not  polydisperse 
particles  for  the  case  rk  =  rt  =  0A  pm.  Fig.  6  shows  that  when  rj 
decreases,  smaller  f-particles  surround  k-  and  1-particles  isolat¬ 
ing  them,  i.e.  breaking  k-1,  k-k  and  1-1  contacts,  resulting  in  lower 
XVTPB  and  narrower  ranges  of  composition  where  XvTpB  >  0.  In  this 
way,  Fig.  6  suggests  to  use  big  pore  former  particles  to  increase 
the  TPB  length.  Nevertheless,  we  must  remember  the  limits  of  our 
model  as  described  in  Section  2.2:  in  order  to  apply  the  percolation 
model,  the  structure  shall  not  collapse  during  the  sintering,  situ¬ 
ation  guaranteed  by  choosing  pore  former  particles  with  proper 
size,  i.e.  similar  to  that  of  k-  and  1-particles.  Thus,  extremely  fine  or 
coarse  pore  former  particles  should  be  avoided.  Instead,  pore  for- 


Volume  fraction  of  k-particles,  y/kas 

Fig.  5.  Mean  hydraulic  radius  as  a  function  of  k-particle  volumetric  composition  for 
different  final  porosities  obtained  by  using  pore  formers  (rfc  =  0.1  |am,  r;  =  0.15  |jim, 
rf=0A  |jim). 
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Fig.  6.  Effect  of  the  radius  of  pore  formers  on  the  TPB  length  per  unit  vol¬ 
ume  obtained  after  sintering  for  a  monosized  mixture  (rfe  =  r/  =  0.1  p,m,  (f^n  =  0.45, 
0  =  15°):  smaller  pore  former  particles  surround  and  isolate  other  particles  resulting 
in  a  decrease  of  TPB  length. 

mers  dimension  should  be  close  to  the  mean  dimension  of  k-  and 
1-particles  or  a  bit  bigger. 

3.4.  Polydisperse  powders 

Concerning  polydisperse  mixtures,  in  the  case  where  pore  for¬ 
mers  are  not  used  (i.e.  we  expect  final  porosity  of  about  0.36)  with 
0=  15°,  we  compare  the  TPB  length  per  unit  volume  as  a  function 
of  composition  after  sintering  among  monosized  and  polydisperse 
mixtures. 

To  describe  polydisperse  powders,  a  symmetrical  numerical 
dispersion  is  considered  as  follows  (described  for  k-particles,  for 
1-particles  it  is  the  same): 

r  rk\  =  ft  “  °*4  •  ft  ftl  =  0-1  •  ft 

rk2  =  ft  ~  0-2  •  ft  Zk2  =  °-2  •  ft 

<  rk3  =  ft  ft3  =  0.4  •  ft  (32) 

rk4  —  rk  +  0.2  •  rk  ft 4  =  0.2  •  ft 

yrk5  =  rk  +  0.4rk  ft5  =  0.1  ■  ft 

while  a  monosized  powder  is  merely  represented  by  its  numerical 
fraction  (e.g.  ft)  and  radius  (e.g.  rk). 

Fig.  7  shows  the  effect  of  the  numerical  dispersion  on  XvTpB  for 
rk  =  ?i  =  0. 1  |jim.  When  a  phase  is  polydisperse  it  behaves  as  it  was 
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Fig.  8.  TPB  length  per  unit  volume  predicted  for  monosized  and  polydisperse  mix¬ 
tures:  symmetrical  volumetric  distribution  =  0.36,  6  =  15°,  rk  =  ?i  =  0.1  pm). 

made  of  bigger  particles,  since  in  a  symmetrical  numerical  distri¬ 
bution  the  major  contribution  in  the  determination  of  numbers  of 
contacts  and  probabilities  of  connection  is  due  to  bigger  particles. 
Fig.  7  can  also  be  read  in  a  different  way:  if  polydisperse  powders 
were  considered  as  monosized  powders,  the  TPB  length  would  be 
overestimated  of  about  25%.  Similar  results  and  considerations  are 
also  reported  by  Kenney  et  al.  [1 1  ],  who  used  computer  simulations 
to  evaluate  XVTPB. 

In  Fig.  8  a  symmetrical  volumetric  distribution,  equal  to  the 
numerical  distribution  described  in  system  (32),  is  used  (i.e.  numer¬ 
ical  fractions  ft/  and  ft  are  replaced  with  volume  fractions  iftj 
and  \j/k).  Fig.  8  shows  little  differences  between  monosized  and 
polydisperse  mixtures  because  the  effects  of  bigger  particles  of  the 
distribution  are  balanced  by  smaller  particles  that  have  the  same 
volume  fraction.  In  this  situation  the  errors  that  we  will  make  on 
the  estimation  of  \PPB  considering  monosized  particles  instead  of 
polydisperse  mixtures  are  low,  but  this  is  principally  due  to  the  fact 
that  the  volumetric  dispersion  is  symmetrical.  Thus,  for  a  generic 
volumetric  dispersion  (i.e.  not  symmetrical)  it  is  reasonable  that 
the  estimation  of  XVTPB  will  be  more  accurate  if  we  consider  that 
distribution  instead  of  using  a  mean  particle  diameter  for  each 
phase. 

4.  Conclusions 
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Fig.  7.  TPB  length  per  unit  volume  predicted  for  monosized  and  polydisperse  mix¬ 
tures:  symmetrical  numerical  distribution  (0^n  =  0.36,  0=  15°,  rk  =  rt  =  0.1  p,m). 


A  percolation  model  has  been  developed  to  predict  effective 
properties  of  composite  electrodes  for  SOFCs  (such  as  TPB  length 
per  unit  volume  and  mean  hydraulic  radius)  from  controllable 
powder  characteristics  such  as  particle  size  distributions  and  com¬ 
positions.  Such  effective  properties  are  needed  to  assist  the  design 
and  optimization  of  electrode  structures  by  continuum-based 
mathematical  models.  The  proposed  model  extends  percolation 
theory  to  multicomponent  polydisperse  mixtures.  As  a  conse¬ 
quence  the  actual  granulometric  distribution  of  powders  may  be 
implemented  in  the  model.  Simulations  show  that  morphological 
properties  of  the  resulting  packing  depend  on  the  granulometric 
distribution  of  the  powders.  Predicted  results  may  vary  consider¬ 
ably  if  the  actual  granulometric  distribution  of  the  powders  is  used 
instead  of  the  mean  particle  dimension,  in  particular  for  asymmet¬ 
rical  or  wide  particle  distributions,  suggesting  that  polydispersion 
of  particles  must  be  considered  in  SOFC  electrode  models. 

Another  feature  introduced  by  the  model  is  the  extension  of  the 
theory  to  account  for  pore  former  particles  (i.e.  porosity),  which 
enables  the  estimation  of  functional  properties  by  the  knowledge 
of  the  amount  of  pore  formers  introduced  before  sintering.  Effective 
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properties  are  strongly  influenced  by  the  variation  of  volume  frac¬ 
tion  of  pore  formers,  which  affects  both  numbers  of  contacts  (i.e. 
probabilities  of  connection)  and  number  of  particles  per  unit  vol¬ 
ume  after  sintering.  The  first  effect  is  generally  neglected  in  other 
percolation  theories  resulting  in  particular  in  an  overestimation 
of  TPB  length.  This  approach  may  be  implemented  for  systems  in 
which  the  cell  design  must  be  expressly  studied  to  keep  the  mass 
transfer  limitations  to  a  minimum  [30]. 

The  model  can  then  be  considered  a  valuable  tool  for  the 
structural  design  of  electrodes,  since  it  provides  a  quantitative 
connection  between  measurable  and  effective  functional  charac¬ 
teristics. 
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Glossary 

avk:  surface  area  exposed  by  k-particles  to  gas  phase  per  unit  volume  (m-1 ) 
avkl:  contact  area  between  connected  k-  and  1-particles  per  unit  volume  (m-1 ) 
av{:  surface  area  exposed  by  1-particles  to  gas  phase  per  unit  volume  (m_  1 ) 
mh:  number  of  particle  sizes  for  a  type  of  particles  ( h  =  k,  1  or  f) 
nv:  number  of  particles  per  unit  volume  (voids  included)  (irr3) 
p:  probability  of  connection  of  a  type  of  particles 
P:  ratio  of  radii  of  particles 

Pmax:  maximum  ratio  of  radii  to  avoid  segregation 
Pmin:  minimum  ratio  of  radii  to  avoid  segregation 
r:  radius  of  a  type  of  particles  (m) 

Rf,:  mean  hydraulic  radius  (m) 

5;  surface  area  fraction  of  a  type  of  particles 
V:  volume  of  a  particle  (m3) 

Zhfh:  overall  average  number  of  contacts  for  all  h-particles  (h  =  k,  1  or  f) 

Zhi:  average  coordination  number  of  hi-particles  ( h  =  k,  1  or  f) 

Zhiy  average  number  of  contacts  between  a  hi-particle  and  f-particles 
Zhig:  average  number  of  contacts  between  a  hi-particle  and  g-particles 
Zhifgj:  average  number  of  contacts  between  a  hi-particle  and  gj-particles 
Z^gj  :  limit  average  number  of  contacts  of  a  hi-particle  with  gj-particles 

Zhi'k  '■  average  number  of  contacts  between  a  hi-particle  and  k-particles 
Zhiy  average  number  of  contacts  between  a  hi-particle  and  1-particles 
Zk:  average  coordination  number  of  k-particles 
Zky.  overall  average  number  of  contacts  for  all  k-particles 
Z\.  average  coordination  number  of  1-particles 
Z\y  overall  average  number  of  contacts  for  all  1-particles 
numerical  fraction  of  a  type  of  particles 
0:  mean  contact  angle  between  two  particles  (rad) 

X”pB:  three-phase  boundary  length  per  unit  volume  (m~2) 

(pfg:  fraction  of  porosity  created  by  pore  formers 
</^n:  final  porosity  (after  sintering) 

0”1;  initial  porosity  (before  sintering) 

Jr:  volume  fraction  of  a  type  of  particles  relative  to  the  total  solid 

Superscripts 

as:  after  sintering 

bs:  before  sintering 

c;  circle  of  contact 

eff:  effective 

kk:  k-k  contact 

Id:  k-1  contact 

Ik:  1-k  contact 

ll:  1-1  contact 

lost:  lost  in  contacts 

overlap:  overlapping  effects  considered 

Subscripts 

f:  pore  former  particles 

fi:  i-th  pore  former  particle  size 

g:  g-particles  (g=  k,  1  or  f  if  not  otherwise  specified) 

gj:  j-th  g-particle  size  (g=  k,  1  or  f  if  not  otherwise  specified) 

h:  h-particles  ( h  =  k,  1  or  f  if  not  otherwise  specified) 

hi:  i-th  h-particle  size  ( h  =  k,  1  or  f  if  not  otherwise  specified) 

k:  k-particles 

ki:  i-th  k-particle  size 

kk:  k-k  contact 

Id:  k-1  contact 

l:  1-particles 

li:  i-th  1-particle  size 

Ik:  1-k  contact 

ll:  1-1  contact 


